{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eddy detection : Antartic Circumpolar Current\n\nThis script detect eddies on the ADT field, and compute u,v with the method add_uv (use it only if the Equator is avoided)\n\nTwo detections are provided : with a filtered ADT and without filtering\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib import style\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\n\npos_cb = [0.1, 0.52, 0.83, 0.015]\npos_cb2 = [0.1, 0.07, 0.4, 0.015]\n\n\ndef quad_axes(title):\n    style.use(\"default\")\n    fig = plt.figure(figsize=(13, 10))\n    fig.suptitle(title, weight=\"bold\", fontsize=14)\n    axes = list()\n\n    ax_pos = dict(\n        topleft=[0.1, 0.54, 0.4, 0.38],\n        topright=[0.53, 0.54, 0.4, 0.38],\n        botleft=[0.1, 0.09, 0.4, 0.38],\n        botright=[0.53, 0.09, 0.4, 0.38],\n    )\n\n    for key, position in ax_pos.items():\n        ax = fig.add_axes(position)\n        ax.set_xlim(5, 45), ax.set_ylim(-60, -37)\n        ax.set_aspect(\"equal\"), ax.grid(True)\n        axes.append(ax)\n        if \"right\" in key:\n            ax.set_yticklabels(\"\")\n    return fig, axes\n\n\ndef set_fancy_labels(fig, ticklabelsize=14, labelsize=14, labelweight=\"semibold\"):\n    for ax in fig.get_axes():\n        ax.grid()\n        ax.grid(which=\"major\", linestyle=\"-\", linewidth=\"0.5\", color=\"black\")\n        if ax.get_ylabel() != \"\":\n            ax.set_ylabel(ax.get_ylabel(), fontsize=labelsize, fontweight=labelweight)\n        if ax.get_xlabel() != \"\":\n            ax.set_xlabel(ax.get_xlabel(), fontsize=labelsize, fontweight=labelweight)\n        if ax.get_title() != \"\":\n            ax.set_title(ax.get_title(), fontsize=labelsize, fontweight=labelweight)\n        ax.tick_params(labelsize=ticklabelsize)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load Input grid, ADT is used to detect eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "margin = 30\n\nkw_data = dict(\n    filename=data.get_demo_path(\"nrt_global_allsat_phy_l4_20190223_20190226.nc\"),\n    x_name=\"longitude\",\n    y_name=\"latitude\",\n    # Manual area subset\n    indexs=dict(\n        latitude=slice(100 - margin, 220 + margin),\n        longitude=slice(0, 230 + margin),\n    ),\n)\ng_raw = RegularGridDataset(**kw_data)\ng_raw.add_uv(\"adt\")\ng = RegularGridDataset(**kw_data)\ng.copy(\"adt\", \"adt_low\")\ng.bessel_high_filter(\"adt\", 700)\ng.bessel_low_filter(\"adt_low\", 700)\ng.add_uv(\"adt\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Identification\nRun the identification step with slices of 2 mm\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "date = datetime(2016, 5, 15)\nkw_ident = dict(\n    date=date, step=0.002, shape_error=70, sampling=30, uname=\"u\", vname=\"v\"\n)\na, c = g.eddy_identification(\"adt\", **kw_ident)\na_, c_ = g_raw.eddy_identification(\"adt\", **kw_ident)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Figures\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kw_adt = dict(vmin=-1.5, vmax=1.5, cmap=plt.get_cmap(\"RdBu_r\", 30))\nfig, axs = quad_axes(\"General properties field\")\ng_raw.display(axs[0], \"adt\", **kw_adt)\naxs[0].set_title(\"Total ADT (m)\")\nm = g.display(axs[1], \"adt_low\", **kw_adt)\naxs[1].set_title(\"ADT (m) large scale, cutoff at 700 km\")\nm2 = g.display(axs[2], \"adt\", cmap=plt.get_cmap(\"RdBu_r\", 20), vmin=-0.5, vmax=0.5)\naxs[2].set_title(\"ADT (m) high-pass filtered, a cutoff at 700 km\")\ncb = plt.colorbar(m, cax=axs[0].figure.add_axes(pos_cb), orientation=\"horizontal\")\ncb.set_label(\"ADT (m)\", labelpad=0)\ncb2 = plt.colorbar(m2, cax=axs[2].figure.add_axes(pos_cb2), orientation=\"horizontal\")\ncb2.set_label(\"ADT (m)\", labelpad=0)\nset_fancy_labels(fig)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The large-scale North-South gradient is removed by the filtering step.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = quad_axes(\"\")\naxs[0].set_title(\"Without filter\")\naxs[0].set_ylabel(\"Contours used in eddies\")\naxs[1].set_title(\"With filter\")\naxs[2].set_ylabel(\"Closed contours but not used\")\ng_raw.contours.display(axs[0], lw=0.5, only_used=True)\ng.contours.display(axs[1], lw=0.5, only_used=True)\ng_raw.contours.display(axs[2], lw=0.5, only_unused=True)\ng.contours.display(axs[3], lw=0.5, only_unused=True)\nset_fancy_labels(fig)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Removing the large-scale North-South gradient reveals closed contours in the\nSouth-Western corner of the ewample region.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kw = dict(ref=-10, linewidth=0.75)\nkw_a = dict(color=\"r\", label=\"Anticyclonic ({nb_obs} eddies)\")\nkw_c = dict(color=\"b\", label=\"Cyclonic ({nb_obs} eddies)\")\nkw_filled = dict(vmin=0, vmax=100, cmap=\"Spectral_r\", lut=20, intern=True, factor=100)\nfig, axs = quad_axes(\"Comparison between two detections\")\n# Match with intern/inner contour\ni_a, j_a, s_a = a_.match(a, intern=True, cmin=0.15)\ni_c, j_c, s_c = c_.match(c, intern=True, cmin=0.15)\n\na_.index(i_a).filled(axs[0], s_a, **kw_filled)\na.index(j_a).filled(axs[1], s_a, **kw_filled)\nc_.index(i_c).filled(axs[0], s_c, **kw_filled)\nm = c.index(j_c).filled(axs[1], s_c, **kw_filled)\n\ncb = plt.colorbar(m, cax=axs[0].figure.add_axes(pos_cb), orientation=\"horizontal\")\ncb.set_label(\"Similarity index (%)\", labelpad=-5)\na_.display(axs[0], **kw, **kw_a), c_.display(axs[0], **kw, **kw_c)\na.display(axs[1], **kw, **kw_a), c.display(axs[1], **kw, **kw_c)\n\naxs[0].set_title(\"Without filter\")\naxs[0].set_ylabel(\"Detection\")\naxs[1].set_title(\"With filter\")\naxs[2].set_ylabel(\"Contours' rejection criteria\")\n\ng_raw.contours.display(axs[2], lw=0.5, only_unused=True, display_criterion=True)\ng.contours.display(axs[3], lw=0.5, only_unused=True, display_criterion=True)\n\nfor ax in axs:\n    ax.legend()\n\nset_fancy_labels(fig)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Very similar eddies have Similarity Indexes >= 40%\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Criteria for rejecting a contour :\n 0. Accepted (green)\n 1. Rejection for shape error (red)\n 2. Masked value within contour (blue)\n 3. Under or over the pixel limit bounds (black)\n 4. Amplitude criterion (yellow)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_a, j_a = i_a[s_a >= 0.4], j_a[s_a >= 0.4]\ni_c, j_c = i_c[s_c >= 0.4], j_c[s_c >= 0.4]\nfig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot (A : {i_a.shape[0]}, C : {i_c.shape[0]} matches)\")\n\nfor i, (label, field, factor, stop) in enumerate(\n    (\n        (\"Speed radius (km)\", \"radius_s\", 0.001, 120),\n        (\"Effective radius (km)\", \"radius_e\", 0.001, 120),\n        (\"Amplitude (cm)\", \"amplitude\", 100, 25),\n        (\"Speed max (cm/s)\", \"speed_average\", 100, 25),\n    )\n):\n    ax = fig.add_subplot(2, 2, i + 1, title=label)\n    ax.set_xlabel(\"Without filter\")\n    ax.set_ylabel(\"With filter\")\n\n    ax.plot(\n        a_[field][i_a] * factor,\n        a[field][j_a] * factor,\n        \"r.\",\n        label=\"Anticyclonic\",\n    )\n    ax.plot(\n        c_[field][i_c] * factor,\n        c[field][j_c] * factor,\n        \"b.\",\n        label=\"Cyclonic\",\n    )\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.plot((0, 1000), (0, 1000), \"g\")\n    ax.set_xlim(0, stop), ax.set_ylim(0, stop)\n    ax.legend()\n\nset_fancy_labels(fig)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}